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ABSTRACT 

It is often assumed that local sources of ionizing radiation have little impact on the 
distribution of neutral hydrogen in the post-reionization Universe. While this is a good 
assumption for the intergalactic medium, analytic arguments suggest that local sources 
may typically be more important than the meta-galactic background radiation for high 
column density absorbers (iVm > 10 17 cm~ 2 ). We post-process cosmological, hydro- 
dynamical simulations with accurate radiation transport to investigate the impact of 
local stellar sources on the column density distribution function of neutral hydrogen. 
We demonstrate that the limited numerical resolution and the simplified treatment of 
the interstellar medium (ISM) that are typical of the current generation of cosmolog- 
ical simulations provide significant challenges, but that many of the problems can be 
overcome by taking two steps. First, using ISM particles rather than stellar particles 
as sources results in a much better sampling of the source distribution, effectively 
mimicking higher-resolution simulations. Second, by rescaling the source luminosities 
so that the amount of radiation escaping into the intergalactic medium agrees with 
that required to produce the observed background radiation, many of the results be- 
come insensitive to errors in the predicted fraction of the radiation that escapes the 
immediate vicinity of the sources. By adopting this strategy and by drastically varying 
the assumptions about the structure of the unresolved ISM, we conclude that we can 
robustly estimate the effect of local sources for column densities A^hi <C 10 21 cm~ 2 . 
However, neither the escape fraction of ionizing radiation nor the effect of local sources 
on the abundance of A^hi > 10 21 cm -2 systems can be predicted with confidence. We 
find that local stellar radiation is unimportant for iVni <C 10 17 cm~ 2 , but that it can 
affect Lyman Limit and Damped Lya systems. For 10 18 < iVni < 10 21 cm~ 2 the im- 
pact of local sources increases with redshift. At z = 5 the abundance of absorbers 
is substantially reduced for Nui 3> 10 17 cm -2 , but at z = the effect only becomes 
significant for Nbi > 10 21 cm -2 . 

Key words: radiative transfer - methods: numerical - galaxies: evolution - galaxies: 
formation - galaxies: high-redshift - intergalactic medium 



1 INTRODUCTION 

After the reionization of the Universe at z > 6, hydrogen re- 
siding in the intergalactic medium (IGM) is kept highly ion- 
ized primarily by the meta-galactic UV background (UVB). 
The UVB is the integrated radiation that has been able to 
escape from sources into the IGM. Because the mean free 
path of ionizing photons is large compared to the scale at 
which the sources of ionizing radiation cluster, the UVB is 
expected to be close to uniform in the IGM. However, close 
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to galaxies, the radiation field is dominated by local sources 
and hence more inhomogeneous. 

Observations of neutral hydrogen (HI) in the Lya for- 
est mostly probe the low-density IGM which is typically far 
from star-forming regions. The statistical properties of the 
Lya forest are therefore in sensitive to the small-s cale fluc- 
tuations in the UVB (e.g., IZuolll992l ; ICroftll2004r i. On the 
other hand, the neutral hydrogen in Damped Lya (DLA; 
i.e., iVm > 10 20 ' 3 cm -2 ) and Lyman Limit systems (LL; i.e., 
10 17 ' 2 < N m < 10 20,3 cm' 2 ), which are thought to originate 
inside or close to galaxies, might be substantially affected by 
radiation fro m local sourc es that are stronger than the am- 
bient UVB (|Gnedinll2O10l) . As a result, the abundances of 
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the high HI column densities may also change significantly 

by locally pr oduced ra diatio n. 

Indeed, ISchavd (2006) and iMiralda-Escudi (|2005l ) 
have used analytic arguments to show that the impact of 
local radiation may be substantial for LL and (sub-)DLA 
systems, but should generally be very small at lower column 
densities. However, relatively little has been done to go 
beyond idealized analytic arguments and to simulate the 
effect of local radiation by taking into account the inho- 
mogeneous distribution of sources and gas in and around 
galaxies. One of the main reasons for this is the compu- 
tational expense of radiative transfer (RT) calculations in 
simulations with large numbers of sources. In addition, high 
resolution is required to capture the distribution of gas on 
small scales accurately. Because of these difficulties, most 
simulations of the cosmological HI distribution have ignored 
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have taken into account local stellar radiatio n but their 
results are inconclu sive: while iNaeamine et alJ (|2010h and 
lYaiima et ail (|2012l ) found that local stellar radiation has a 
neglig ible impact on the distribution of HI, Fumagall i et alJ 
l)201ll ) found that the HI column density distribution above 
the Lyman limit is reduced by ~ 0.5 dex due to local stellar 
radiation. 

In this paper, we investigate the impact of local stellar 
radiation on the HI distribution by combining cosmological 
hydrodynamical simulations with accurate RT. We find that 
the inclusion of local sources can dramatically change the 
predicted abundance of strong DLA systems and, depend- 
ing on redshift, also those of LL and weak DLA systems. 
On the other hand, lower HI column densities are hardly 
affected. We also show that resolution effects have a major 
impact. For instance, the resolution accessible to current cos- 
mological simulations is insufficient to resolve the interstel- 
lar medium (ISM) on the scales relevant for the propagation 
of ionizing photons. The limited resolution also affects the 
source distribution which can change the resulting HI dis- 
tribution, especially in low-mass galaxies. On top of that, 
assumptions about the structure of the unresolved multi- 
phase ISM significantly affect the escape of stellar radiation 
into the IGM. Therefore, any attempt to use cosmological 
simulations to investigate the impact of local stellar radia- 
tion on the HI distribution may suffer from serious numerical 
artifacts. 

Some of these difficulties can be circumvented by tun- 
ing the luminosities of the sources such that the escaped 
radiation can account for the observed UVB. Then the in- 
teraction between the radiation that reaches the IGM and 
the intervening gas can be studied on scales that are prop- 
erly resolved in the simulations. We adopt this procedure to 
generate the observed UVB for various ISM models but find 
that our fiducial simulation reproduces the observed UVB 
without any tuning. 

Among the known sources of radiation, quasars and 
massive stars are the most efficient producers of hydro- 
gen ionizing photons and are therefore thought to be the 
main contributors to the UVB. Star- forming galaxies how- 
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iFaucher-Giguere et al.ll2008ah . Therefore, we only account 
for local radiation that is produced by star formation in our 
simulations. Moreover, we assume that the ionizing emissiv- 
ity of baryons strictly follows the star formation rate. We use 
star-forming gas particles rather than young stellar particles 
as ionizing sources. Since the gas consumption time scale in 
the ISM is much larger than the lifetime of massive stars 
(~ 10 9 yr vs. ~ 10 r yr), there are many more star- forming 
gas particles than there are stellar particles young enough 
to efficiently emit ionizing radiation. Therefore, using star- 
forming gas particles as ionizing sources allows us to sample 
the source distribution better and hence to reduce the im- 
pact of the limited resolution of cosmological simulations. 

The structure of the paper is as follows, in iJH we use 
the observed relation between the star formation rate and 
gas surface densities to provide an analytic estimate for the 
photoionization rate that is produced by young stars in a 
uniform and optically thick ISM that is in good agreement 
with our simulation results. In 331 we discuss the details of 
our numerical simulations and RT calculations. The simu- 
lation results are presented in $4] We show how the local 
stellar radiation can generate the observed UVB. We also 
discuss the effect of unresolved ISM on the escape fraction 
of ionizing radiation and we investigate the impact of lo- 
cal stellar radiation on the HI column density distribution. 
Finally, we conclude in ij5] 



2 PHOTOIONIZATION RATE IN 
STAR-FORMING REGIONS 

A strong correlation between gas surface density and star 
form ation rate has been observed in lo w-redshift galaxies 
(e.g., lKennicutt|[l998l ; iBigiel et al.l [20081 ). In principle, such 
a relation can be combined with simplified assumptions to 
derive a typical photoionization rate that is expected from 
stellar radiation. In this section, we use this approach to 
estimate the average ionization rate that is produced by star 
formation as a function of gas (surface) density on galactic 

scales. 

Assuming a IChabrierl ([2003) initial mass function 
(IMF), the observed Kennicutt-Schmidt law provides a re- 
lation between gas surface density, E gas , and star formation 
rate surface density, £*, on kilo-parsec scales ( Kennicuttl 
1 19981) : 



E* « 1.5 x 10~ 4 M yr _1 kpc" 



lM pc- 



(1) 



Equation |T]) can be used to derive a relation between ioniz- 
ing emissivity (per unit area), E T , and the gas surface den- 
sity. As we will discuss in 33.41 stellar population synthesis 
models indicate that the typical number of ionizing photons 
produced per unit time by a constant star formation rate is 
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Furthermore, as both models and observations suggest, the 
escape fraction of ionizing photons from galaxies is <C 1 (e.g. , 
Shaplev et al.ll2006l;lGnedin et al.ll2008l;IVanzella et alfcOlfj; 
Yaiima et al.l 20121 ; iPaardekooper et al.l l201ll ~ Kim et al.1 
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12012ft . This allows us to assume that most of the ionizing 
photons that are produced by star-forming gas are absorbed 
on scales < kpc. Therefore, the hydrogen photoionization 
rate, on galactic scales, can be computed: 



(3) 



where JVh is the hydrogen column density which can be ob- 
tained from the gas surface density: 



N H w 9.4 x 10 19 cm~ 2 
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where X is the hydrogen mass fraction. After assuming X = 
0.75 and substituting equations {TJ, ([2]) and Q in equation 
((3]), one gets 



8.5 x 10~ 14 s _1 



10 21 cm- 
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If the scale height of the disk is similar to the local Jeans 
scale, the hydrogen column density can be computed as a 
function of the hydrogen numb er density (|Schavel200ll , l2004l : 
ISchave fc Dalla Vecchiall200Sl ) 



N H ~ 2.8 x 10 21 cm -2 
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where T4 = T/10 4 K, / g is the total mass fraction in gas 
and fth is the fraction of the pressure that is thermal (i.e., 
Pth = /th-Prcr). In addition, for deriving equation ((6| we 
have adopted an adiabatic index, 7 = 5/3, mean particle 
mass, jj, = 1.23 nm. After eliminating A^h between equation 
([5]) and equation ([6]), the photoionization rate can be written 
as a function of the hydrogen number density 
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Based on equation l[7[l. the photoionization rate pro- 
duced by star formation is only weakly sensitive to the tem- 
perature, total gas fraction, / g , and the fraction of the pres- 
sure that is thermal, / t h- The photoionization rate in equa- 
tion Q is also very weakly dependent on the gas density. As 
we will discuss in N4.ll our simulations show the same trend 
and are in excellent agreement with this analytic estimate. 
However, one should note that due to the simplified assump- 
tions we adopted to derive equation (0), it provides only an 
order-of-magnitude estimate for the typical photoionization 
rate that is produced by young stars in the ISM, on kilo- 
parsec scales. This relation does not capture the inhomo- 
geneity of the ISM that may result in large fluctuations in 
the radiation field on smaller scales. We note that the rela- 
tion we found here between the gas density and stellar pho- 
toionization rate is a direct consequence of the underlying 
star formation law and is independent of the total amount 
of star formation in a galaxy. 



3 SIMULATION TECHNIQUES 

In this section we describe different parts of our simulations. 
We start with discussing the details of the hydrodynamical 
simulations that are post-processed with RT calculations. 
Then we explain our RT that includes the radiation from 



local stellar radiation, the UVB and recombination radia- 
tion. 



3.1 Hydrodynamical simulations 

We use a set of cosmological simulations that are performed 
using a modified and extended version of the smoothed par- 
ti cle hydrodynam ics (SPH) code GADGET-3 (last described 
in lSpringelll2005l V The physical processes that are included 
in the simulations are identical to what has been used in the 
reference simulation of the Overwhelmingly L arge Simula- 
tions (OWLS) described in lSchave etafl (jboicT) . Briefly, we 
use a s ubgrid pressure-dependen t star formation prescrip- 
tion of ISchave fc Dalla Vecchial (| 200 8) which reproduces 
the observed Kenni cutt-Schmidt law. We use the chemo- 
dynamics model of W lcrsma et all l|2009bh which follows 
the a bundances of eleven elements assuming a IChabrierl 
(2003) IMF. These abundances are used for calculating ra- 
diative cooling/heating rates, element-by-element and in the 
pre sence of the uniform cosm ic microwave background and 
the lHaardt fc Madaul (I200ll1 UVB model dWiersma et ail 
l2009al ). We note that the lHaardt fc Madaul |200ll ) UVB 
model has been shown to be consistent with observations 
of HI (|Altav et all l201ll: IRahmati et all 120121 ') and metal 



absorption lines ( Aguirre et al. 20081 ') . We model galactic 
winds driven by star formation using a kinetic feedback 
recipe that assumes 40% of the kinetic energy generated 
by Type II SNe is injected as outflows with initial velocity 
of 600 kms~ and with a ma ss loading parameter r] — 2 
jDalla Vecchia fc Schavell200St ). 

We adopt cosmological parameters consistent with the 
most recent WMAP year 7 results: {fi m = 0.272, tt h = 
0.0455, fl A = 0.728 a 8 = 0.81, n s = 0.967, h = 0.704} 
l|Komatsu et al.ll201ll ). Our reference simulation has a peri- 
odic box of L = 6.25 comoving /i -1 Mpc and contains 128 3 
dark matter particles with mass 6.3 x 10 6 /i _1 M Q and an 
equal number of baryons with initial mass 1.4 x 10 6 /i _1 M . 
The Plummer equivalent gravitational softening length is set 
to e com = 1.95 /i _1 kpc and is limited to a minimum physical 
scale of e P rop = 0.5 ft _1 kpc. We also use a simulation with a 
box size identical to our reference simulation but with 8 (2) 
times better mass (spatial) resolution to assess the effect of 
resolution on our findings (see Appendix IB2|) . 

In our hydrodynamical simulations, ISM gas particles 
(which all have densities n H > 0.1 cm -3 ) follow a poly tropic 
equation of state that defines their temperatures. These tem- 
peratures are not physical and only me asure the imposed 
pressure (jSchave fc Dalla Vec chia 200$). Therefore, when 
calculating recombination and collisional ionization rates, 
we set the temperature of ISM particles to Tism = 10 4 K 
which is the typical temperature of the warm-neutral phase 
of the ISM. Furthermore, we simplify our RT calculations 
by assuming that helium and hydrogen absorb the same 
amount of ionizing radiation per unit mass and by ignoring 
dust absorption a nd the possibility tha t some hydrogen may 
be molecular (see lRahmati et aLll2012j for more discussion). 



3.2 Radiative transfer 

For the RT c alcula tio ns w e use TRAPH IC (se e 
Pawli k fc Schavd 120081 . l201ll : IRahmati etail l2012h . 



© 0000 RAS, MNRAS 000, 000-000 



4 A. Rahmati et al. 



TRAPHIC is an explicitly photon-conserving RT method 
designed to exploit the full spatial resolution of SPH sim- 
ulations by transporting radiation directly on the irregular 
distribution of SPH particles. 

The RT calculation in TRAPHIC starts with source par- 
ticles emitting photon packets to their neighbors. This is 
done using a set of JVec tessellating emission cones, each 
subtending a solid angle of Ah/N-ec- The propagation direc- 
tions of the photon packets are initially parallel to the cen- 
tral axes of the emission cones. In order to improve the an- 
gular sampling of the RT, the orientations of these emission 
cones are randomly rotated between emission time steps. 
We adopt Aec = 8 for computational efficiency. However, 
we note that our results are insensitive to the precise value 
of A^ec, thanks to the random rotations. 

After emission, photon packets travel along their propa- 
gation direction from one SPH particle to its neighbors. Only 
neighbors that are inside transmission cones can receive pho- 
tons. Transmission cones are defined as regular cones with 
opening solid angle An /Ntc, and are centered on the prop- 
agation direction. The parameter Ntc sets the angular res- 
olution of the RT and we adopt Ntc = 64 which produce s 
converged results (see Appendix CI of lRahmati et al.ll2012h . 
To guarantee the independence of the RT from the distri- 
bution of SPH particles, additional virtual particles (ViP) 
are introduced. This is done wherever a transmission cone 
contains no neighboring SPH particle. ViPs do not affect the 
underlying SPH simulation and are deleted after the photon 
packets are transferred. 

Furthermore, arriving photon packets are merged into 
a discrete number of reception cones, ATrc = 8. This makes 
the computational cost of RT calculations with TRAPHIC in- 
dependent of the number of radiation sources. This feature is 
particularly important for the purpose of the present work as 
it enables the RT calculations in cosmological density fields 
with large numbers of sources. 

Reception cones are also used for emission from gas par- 
ticles (e.g., ionizing radiation from star-forming gas parti- 
cles, recombination radiation). This reduces computational 
expenses while yielding accurate results. Photon packets are 
isotropically emitted into reception cones which are already 
in place at each SPH particle. This obviates the need for 
constructing any additional emission cone tessellations. Us- 
ing this recipe, the emission of radiation by star-forming gas 
particles is identical to that of the emission of diffuse radia- 
tion by recombining gas particles, and is described in more 
detail in Raicevic et al. (in prep.). 

Photon packets emitted by the three main radiation 
sources we consider in our study (i.e., UVB, RR and stellar 
radiation) are channeled in separate frequency bins and are 
not merged with each other. This enables us to compute the 
contribution of each component to the total photoionization 
rate. The total amount of the absorbed radiation, summed 
over all frequency bins, determines the ionization state of 
the absorbing SPH particles. The time-dependent differen- 
tial equations that control the evolution of the ionization 
states of different species (e.g., H, He) are solved using a 
sub-cycling scheme that allows us to choose the RT time 
step independent of the photoionization and recombination 
times. 

In our RT calculations we use a set of numer i cal pa - 
rameters identical to that used in iRahmati etlll (|2012h . 



54 














Z = 


1.6 Z 

© 


1 M vr"' 






z = 


0.66 Z 

© 








z = 


0.33 Z 




1 




z = 






to 
Of 










o 

eo" 
j3 






/- 


Kroupa IMF 
Salpeter IMF 


52 




' 

** 


■ ' 





10 5 10 6 10 7 10' 

Time (yr) 



Figure 1. The number of hydrogen ionizing photons produced for 
a constant star formation rate of 1 Mq yr -1 as a function of time 
since the onset of star formation. These results are calculated us- 
ing STARBURST99. Red, green, blue and purple curves indicate 
metallicities of Z/Z = 1.6, 6.6 x 1CT 1 3.3 x 10" 1 , &3.3xl0~ 2 , 
respectively. Solid and dashed curves show results for the Kroupa 
and Salpeter IMF, respectively. All the curves with different 
metallicities and IMFs converge to similar equilibrium values at 
Q 7 ~ 2 x 10 53 photons per second, after about ~ 5 — 10 Myr. 

These parameters produce converged results. In addition to 
the parameters mentioned above, we use 48 neighbors for 
SPH particles, 5 neighbors for ViPs and an RT time step 
of At = 1 Myr ^YTz)' I° mzm g photons are propagated at 
the speed of light, inside the simulation box with absorbing 
boundaries, until the equilibrium solution for the hydrogen 
neutral fractions is reached. To facilitate this process, the 
RT calculation is starting from an initial neutral state, ex- 
cept for the gas at low densities (i.e., n gas < 1 x 10~ 3 cm -3 ) 
and high temperatures (i.e., T > 10 J K) which is assumed 
to be in equilibrium with the UVB photoionization rate and 
its collisional ionization rate. Typically, the average neutral 
fraction in the simulation box does not evolve after 2-3 light- 
crossing times (the light-crossing time for the Lbox = 6.25 
comoving /i -1 Mpc is « 7.2 Myr at z = 3). 

We continue this section by briefly describing the im- 
plementation of the UVB, diffuse recombination radiation 
and stellar radiation. 



3.3 Ionizing background radiation and diffuse 
recombination radiation 

In principle, local sources of ionization inside the simulation 
box should be able to generate the UVB. However, the box 
size of our simulations is smaller than the mean free path 
of ionizing photons which makes the simulated volume too 
small to generate the observed UVB intensity (see H4.3p . In 
addition, a considerable fraction of the UVB at z < 3 is pro- 
duced by quasars which are not included in our simulations. 
Therefore, we impose an additional UVB in our simulation 
box. 

The implementation of the UVB is identical to that of 
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Figure 2. The HI column density distribution of a Milky Way-like galaxy in the reference simulation at z = 3 (top) and 2 = (bottom). 
The left panels show the total hydrogen column density distribution. The middle panels show the HI column densities in the presence 
of collisional ionization, photoionization from the UVB and recombination radiation. In the right panel, the HI distribution is shown 
after photoionization from local stellar radiation is added to the other sources of ionization. The images are 100 h _1 kpc (proper) on a 
side. Comparison between the right and middle panels shows that local stellar radiation substantially changes the HI column density 



iRahmati et"aL1 l|2012tl . The ionizing background radiation is 
simulated as plain-parallel radiation entering the simulation 
box from its sides and the injection rate of the UVB ionizing 
photons is normalized to the desired photoionization rate in 
the absence of any absorption (i.e., the so-called optically 
thin limit). 

We set the effective photoionization rate and spectral 
shape of the UVB radiation at different re dshifts based on 
the UVB model of lHaardt fc Madaul | 20011 ) for quasars and 
galaxies. The same UVB model has been used to calcu- 
late heating/cooling in our hydrodynamical simulations and 
has been shown to be consistent with observations of HI 
(Altav ct al. 2011; Rahmati ct al. 2012) and metal absorp- 
tion lines 1 Aguirre et alj 2008) at z ~ 3. We adopt the gray 
approximation instead of an explicit multifrequency treat- 
ment for the UVB radiation. Our tests show that a multi- 
frequency treatment of the UVB radiation does not signifi- 
cantly change t he resulting hydrogen neutral fractions (see 
Appendix Dl of lRahmati et al.ll2012T ) 

In addition, hydrogen recombination radiation (RR) is 
simulated by making all SPH particles isotropic radiation 
sources with emissivities based on their recombination rates. 
We do not account for the redshifting of the recombination 



photons and assume that they are monochromatic with en- 
ergy 13.6 eV (see Raicevic et al. in prep.). 



3.4 Stellar ionizing radiation 

The ionizing photon production rate of star-forming galax- 
ies is dominated by young and massive stars which have 
relatively short life times. In cosmological simulations with 
limited mass resolutions, the spatial distribution of newly 
formed stellar particles (e.g., with ages less than a few tens 
of Myr) may sample the locally imposed star formation rates 
relatively poorly. As we show in H4.21 the distribution of star 
formation is better sampled by star-forming SPH particles 
and this is particularly important for low and intermediate 
halo masses. For this reason, we use star-forming particles 
as sources of stellar ionizing radiation. 

For a constant star formation rate, the production rate 
of ionizing photons reaches an equilibrium value within 
~ 5 — 10 Myr (i.e., the typical life-time of massive stars). This 
equilibrium photon production rate per unit star formation 
rate can be calculated using stel lar population synth esis 
models. We used STARBURST 99 (iLeitherer et al.|[l999Tl to 
calculate this emissivity for the lKroupa ( 200ll ) and Salpeterl 



© 0000 RAS, MNRAS 000, 000-000 



6 A. Rahmati et al. 



(1955) IMF and a metallicity consistent with the typical 
metallicities of star-forming particles in our simulations (the 
median metallicity of star-forming particles in our simula- 
tions are within the range 10 _1 < Z/Z Q < 1 at redshifts 
0-5). We found that for a constant star formation rate 
of 1 Mq yr _1 and after ~ 10 Myr, the photon production 
rate of ionizing photons converges to Q 7 ~ 2 x 10 53 s , 
which is the value we used in equation ([2]). As shown in 
Figure [T] the ionizing photon production rate varies only by 
< ±0.2 dex if the metallicity changes by ~ ±1 dex. Also, 
reasonable variations in the IMF (e.g., the Chabrier IMF) 
do not significantly change our adopted value for the pho- 
ton production rate. For the spectral shape of the stellar 
radiation, we adopt a blackbody spectrum with tempera- 
ture Tbb = 5x 10 4 K which is consistent with the spectrum 
of massive young stars. 



4 RESULTS AND DISCUSSION 

In this section, we report our findings based on RT simu- 
lations that include the UVB, recombination radiation and 
local stellar radiation. The RT calculations are performed 
by post-processing a hydrodynamical simulation with 128 3 
SPH particles in a 6.25 comoving /i~ 1 Mpc box and at red- 
shifts z = — 5. In ij4.ll we compare different photoionizing 
components and assess the impact of local stellar radiation 
on the HI distribution. In £14.21 we illustrate the importance 
of using star-forming SPH particles instead of using young 
stellar particles as ionizing sources. In £|4.3I we show that in 
our simulations, the predicted intensity of stellar radiation 
in the IGM is consistent with the observed intensity of the 
UVB and calculate the implied average escape fraction. In 
i)4.4l we show that assumptions about the unresolved prop- 
erties of the ISM are very important and we discuss the 
impact of local stellar radiation on the HI column density 
distribution. 

4.1 The role of local stellar radiation in hydrogen 
ionization 

As we showed in iRahmati et~ai] (|2012T ). the UVB radiation 
and collisional ionization are the dominant sources of ion- 
ization at low densities where the gas is shock heated and 
optically thin. However, self-shielding prevents the UVB ra- 
diation from penetrating high-density gas. In addition, much 
of the ISM has temperatures that are too low for collisional 
ionization to be efficient. Therefore, because UVB photoion- 
ization and collisional ionization are both inefficient in those 
regions, the ionizing radiation from young stars becomes the 
primary source of ionization. 

Figure [2] shows the distribution of neutral hydrogen in 
and around a galaxy in our reference simulation at redshifts 
z = 3 (top row) and (bottom row) and illustrates that the 
distribution of HI in galaxies may change significantly as a 
result of local stellar radiation. The total mass of the halo 
hosting this galaxy at z = 3 and is M200 = 2 x 10 x1 Mq and 
1.1 x 10 12 Mq, respectively. Comparing the panels in the left 
and middle columns, we see that while collisional ionization 
and photoionization by the UVB ionize the low-density gas 
around the galaxy, they do not change the high HI column 
densities in the inner regions. Comparing the yellow regions 



in the right and middle panels of Figure [2] demonstrates that 
local stellar radiation significantly changes the distribution 
of the HI at N m ~ 10 21 cm" 2 . 

In the left panel of Figure [3] the contributions of the 
photoionization rates from different radiation components 
are plotted against the total hydrogen number density at 
z — 3. The purple solid line shows the total photoionization 
rate. The blue dashed, green dotted and red dot-dashed lines 
show respectively the contribution of the UVB, recombina- 
tion radiation (RR) and local stellar radiation (LSR). The 
resulting hydrogen neutral fraction as a function of gas den- 
sity is shown by the purple solid curve in the right panel 
of Figure [3] which illustrates the significant impact of LSR 
at high densities. We note that there is a sharp feature in 
the hydrogen neutral fraction at 10 _1 < n H < 10 -2 cm -3 
in the presence of local stellar radiation. This feature corre- 
sponds to a sharp drop-off in the photoionization rate from 
local stellar radiation at the same densities (see the red dot- 
dashed curve in the left panel of Figure [3]) . The left panel 
of Figure [4] shows photoionization rate profiles in spheri- 
cal shells around haloes with 10.5 < log 10 M200 < 11 Mq at 
z = 3. The right panel of Figure[4]shows the resulting neutral 
hydrogen fraction profile (purple solid curve) which is com- 
pared with the simulation that does not include local stellar 
radiation (green dashed curve). We note that the trends in 
the photoionization rate and hydrogen neutral fraction pro- 
files do not strongly depend on the chosen mass bin as long 
as log 10 M200 > 9.5 M Q . 

As Figures [3] and |4] show, at low densities (i.e., at large 
distances from the centers of the halos) the gas is highly 
ionized by a combination of the UVB and collisional ioniza- 
tion and this optically thin gas does not absorb a significant 
fraction of the ambient ionizing radiation. Consequently, at 
n H < 10~ 4 cm -3 , which corresponds to typical distances 
R > 1 comoving Mpc from the centers of the haloes, the 
UVB photoionization rate does not change with density. 
By increasing the density, or equivalently by decreasing the 
distance to the center of the haloes, the optical depth in- 
creases and eventually the gas becomes self-shielded against 
the UVB. This causes a sharp drop in the UVB photoioniza- 
tion rate at the self-shielding density (see the blue dashed 
curves in the l e ft pan els of Figures [3] and [4]) . As shown in 
IRahmati et"ai1 (|2012h . the UVB photoionization rate shows 
a similar behavior in the absence of local stellar radiation. 
However, local stellar radiation increases the photoioniza- 
tion rate at high and intermediate densities. This decreases 
the hydrogen neutral fractions around the self-shielding den- 
sities, allowing the UVB radiation to penetrate to higher 
densities. Consequently, the local stellar radiation increases 
the effective self-shielding threshold against the UVB radi- 
ation slightly (by ~ 0.1 dex) compared to the simulation 
without the radiation from local stars (not shown). 

As shown by the red dot-dashed curve in the left panel 
of Figure|3] at densities 0.1 < n H < 1 cm~ 3 the median pho- 
toionization rate due to local stellar radiation increases with 
decreasing density. The main reason for this is the superpo- 
sition of radiation from multiple sources (i.e., star-forming 
SPH particles) as the mean free path of ionizing photons 
increases with decreasing density (see Figure IB1|I . This is 
also seen in the left panel of Figure [3] as increasing stellar 
photoionization rate with increasing the distance from the 
center of halos for R < 1 comoving kpc. On the other hand, 
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Figure 3. Left: Photoionization rate as a function of density due to different radiation components in the reference simulation at z = 3. 
The purple solid curve shows the total photoionization rate. The blue dashed, green dotted and red dot-dashed curves show respectively 
the photoionization rates due to the UVB, diffuse recombination radiation (RR) and local stellar radiation (LSR). Right: The hydrogen 
neutral fraction as a function of density for the same simulation is shown with the purple solid line. For comparison, the results for 
a simulation without local stellar radiation (green dashed curve) and a simulation with the UVB radiation in the optically thin limit 
(i.e., no absorption; brown dotted curve) are also shown. The curves show the medians and the shaded areas around them indicate the 
f5% — 85% percentiles. HI column densities corresponding to each density in the presence of all ionization sources are shown along the 
top x-axis. The photoionization due to local stellar radiation exceeds the UVB photoionization rate at high densities and compensates 
the effect of self-shielding. This produces a hydrogen neutral fraction profile that is very similar to what is expected from the UVB in 
the optically thin limit. 
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Figure 4. Left: Median photoionization rate profiles (comoving) due to different radiation components for haloes with 10.5 < 



log 10 M 2 oo < If M at 



3. Blue solid, green dotted and red dot-dashed curves show the photoionization rates of the UVB, re- 



combination radiation (RR) and local stellar radiation (LSR). Right: Median hydrogen neutral fraction profiles for the same halos with 
local stellar radiation (purple solid curve), without local stellar radiation (green dashed curve) and with only the UVB radiation in the 
optically thin limit (dotted brown). In both panels, the vertical dotted line indicates the median i?200 radius of the halos in the chosen 
mass bin. The shaded areas around the medians indicate the 15% — 85% percentiles. The top axis in each panel shows the median 
density at a given comoving distance from the center of the haloes. The photoionization due to local stellar radiation exceeds the UVB 
photoionization rate close to galaxies and compensates the effect of self-shielding. This produces a hydrogen neutral fraction profile that 
is very similar to what is expected from the UVB in the optically thin limit. 
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at densities lower than the star formation threshold (i.e., 
n H < 0.1 cm -3 ), the gas is typically at larger distances from 
the star-forming regions. Therefore, the photoionization rate 
of local stellar radiation drops rapidly with decreasing den- 
sity. The star formation rate density averaged on larger and 
larger scales becomes increasingly more uniform. This causes 
the photoionization from galaxies that are emitting ionizing 
radiation to produce a more uniform photoionization rate at 
the lowest densities. Note that at the highest densities, the 
photoionization rate from local stellar radiation agrees well 
with the analytic estimate presented in iJU 

In the absence of local stellar radiation, recombina- 
tion radiation pho toionization rate pea ks at around the self- 
shielding density |Rahmati et al.lfeoij '). The reason for this 
is that the production rate of ionizing photons by recombi- 
nation radiation is proportional to the density of the ion- 
ized gas. At densities below the self-shielding threshold, hy- 
drogen is highly ionized and recombination rate is propor- 
tional to the total hydrogen density. At much higher den- 
sities, the gas is nearly neutral and little recombination ra- 
diation is generated. As the left panel in Figure [3] shows, 
at high densities the situation changes dramatically if we 
include local stellar radiation. Since the gas at high densi- 
ties (e.g., n H 2> 10 -2 cm -3 ) is optically thick, recombina- 
tion radiation ionizing photons produced at these densities 
are absorbed locally. In equilibrium, recombination radiation 
photoionization rate at densities above the self-shielding is 
therefore a constant fraction of the total ionization rat^S 
Below the self-shielding density, the recombination radia- 
tion photoionization rate decreases with decreasing density 
and asymptotes to a background rate. However, in reality 
recombination photons cannot travel to large cosmological 
distances without being redshifted to frequencies below the 
Lyman edge. Therefore, our neglect of the cosmological red- 
shifting of recombination radiation leads us to overestimate 
the photoionization rate due to recombination radiation on 
large scales. On the other hand, because of the small size 
of our simulation box, the total photoionization rate that 
is produced by recombination radiation remains negligible 
compared to the UVB photoionization rate and neglecting 
the redshifting of recombination radiation is not expected 
to affect our results. 

The purple solid curve in the right panel of Figure 
[3] shows the hydrogen neutral fractions in the presence of 
the UVB, recombination radiation and local stellar radi- 
ation. For comparison, the hydrogen neutral fractions in 
the absence of local stellar radiation, and for the optically 
thin gas that is photoionized only by the UVB, are also 
shown (with the green dashed and brown dotted curves re- 
spectively). Hydrogen at densities n H > 10 -1 cm -3 is self- 
shielded and mostly neutral if the UVB and recombination 
radiation are the only sources of photoionization. However, 
local stellar radiation significantly ionizes the gas at inter- 
mediate and high densities. As mentioned in the typical 

1 This can be explained noting that in equilibrium, recombi- 
nation and ionization rates are equal. Depending on temper- 
ature, Ri 40% of recombination photons are ionizing photons 
which are absorbed on th e spot in optically thick gas (e.g., 
lOsterbrock fc Ferlandll2006h . Therefore, the photoionization rate 
produced by recombination radiation is also ss 40% of the total 
ionization rate in dense and optically thick gas with T ~ 10 4 K. 



photoionization rate that is produced by star-forming gas is 
Tsf ~ 10 -13 s - , which is comparable to the photoioniza- 
tion rate of the unattenuated UVB at z ~ 0. Consequently, 
the hydrogen neutral fractions in the presence of local stellar 
radiation are much closer to the optically thin case. It is also 
worth noting that the scatter around the median hydrogen 
neutral fractions is largest for intermediate densities (i.e., 
10 -3 cm -3 < n H < 1cm -3 . This is closely related to the 
large scatter in the photoionization rate produced by local 
stellar radiation (see the left panel in Figure [3]). At these 
densities, the large scatter in the distances to the nearest 
sources and RT effects like shadowing produce a large range 
of photoionization rates. For the gas at very high and very 
low densities on the other hand, the scatter becomes smaller 
because the relative distribution of sources with respect to 
the absorbing gas becomes more uniform. 

The trends discussed so far are qualitatively similar at 
redshifts other than z — 3: Although the star formation 
rates evolve at high densities, the photoionization rate due 
to local stars is set by the underlying star formation law (see 
SJ2| which does not change with time. The peak of the stel- 
lar photoionization rate at n H ~ 10 -1 cm" 3 , produced by 
the superposition of multiple sources, exists at all redshifts 
but it moves to slightly higher densities at lower redshifts. 
At densities immediately below the adopted star formation 
threshold (i.e., n H < 0.1 cm -3 ), the photoionization rate 
produced by local sources drops rapidly with decreasing den- 
sity. The distribution of star formation in the simulation box 
is almost uniform on large scales. Therefore, at the lowest 
densities, local stellar radiation produces a photoionization 
rate which is not changing strongly with density. However, 
at low redshifts (e.g., z < 1), the star formation density de- 
creases significantly and becomes highly non-uniform on the 
scales probed by our small simulation box. Consequently, 
the resulting photoionization rate due to young stars does 
not converge to a constant value at low densities at these 
redshifts. For the same reason, at low densities the scat- 
ter around the median stellar photoionization rate increases 
with decreasing redshift (not shown). As we will discuss in 
i|4.3l if we correct for the small size of our simulation box, 
the asymptotic photoionization rate due to stellar radiation 
that has reached the IGM, is consistent with the intensity 
of the observed UVB at the same redshift. 

4.2 Star-forming particles versus stellar particles 

In cosmological simulations, the integrated instantaneous 
star formation rate of the simulation box closely matches 
the total amount of mass converted into stellar particles. 
However, due to limited resolution, the spatial distribution 
of young stellar particles (e.g., those with ages ~ 10 Myr) 
may not sample the spatial distribution of star formation in 
individual haloes very well. 

This issue is illustrated in the left panel of Figure [5J 
where the blue crosses show the total instantaneous star for- 
mation rates inside galaxies in our reference simulation at 
z — 3. These star formation rates are calculated by the sum- 
mation of the star formation rates of all SPH particles in a 
given galaxy. The star formation rates can also be estimated 
by measuring the average rate by which stellar particles are 
formed in a given galaxy over some small time interval. The 
star formation rate calculated using this method (i.e., mea- 
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Figure 5. Star-formation activity of haloes in the reference simulation (6.25 Mpc , 128 3 ) at z = 3. In the left panel the blue crosses 
show the total instantaneous star formation rate for a given halo computed using the gas particles and the red circles show the star 
formation rate calculated by dividing the total stellar mass formed during the last 20 Myr by 20 Myr. Haloes with zero star formation 
rates are shown in the bottom of the plot. For massive haloes the two measures agree but as a result of limited mass resolution and the 
stochastic nature of the star formation algorithm, they start to differ substantially for haloes with M200 10 10 Mq . A large fraction of 
low-mass haloes does not contain any young stellar particles and the median star formation rate calculated using young stellar particles 
drops to zero. The right panel shows the number of ionizing sources for different halos in the same simulation. While the blue crosses 
show the number of star-forming gas particles in each halo, red circles indicate the number of young stellar particles (i.e., younger than 
20 Myr). In both panels the horizontal dotted lines correspond to a single stellar particle and the blue dashed and red dotted curves 
show respectively the medians for star-forming gas and young stellar particles. 
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Figure 6. Median photoionization rate profiles (comoving) for local stellar radiation emitted by star-forming SPH particles (SF-SPH; 
blue solid curves) and young (< 20 Myr) stellar particles (Y Stars; red dotted curves) at 2 = 3. The left panel shows the photoionization 
rate profiles for a halo with M200 = 2 X 10 11 Mq while the right panel shows the same profiles for haloes with 10 < log 10 M200 < 10-5 Mq . 
In both panels the vertical dotted line indicates the median -R200 radius of the halos in the illustrated mass bin. The shaded areas around 
the medians indicate the 15% — 85% percentiles. The top axis in each panel shows the median density at a given comoving distance from 
the centers of the haloes. While the photoionization rate profiles produced by star-forming gas particles and young stellar particles are 
similar for massive haloes (left), they are substantially different for haloes with M200 $ 10 10 Mq (right). 
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suring the rate by which stellar particles are formed during 
the last 20 Myr) , is indicted in the left panel of Figure [5] 
by the red circles. For massive galaxies, with high star for- 
mation rates, the formation rate of stellar particles agrees 
reasonably well with the star formation rates computed from 
the gas distribution. However, most haloes with < 10 3 SPH 
particles (i.e., M200 <• 10 10 M ) do not contain any young 
stellar particles despite having non-zero instantaneous star 
formation rates. Consequently, for the few low-mass galaxies 
that by chance contain one or more young stellar particles, 
the implied star formation rates are much higher than the 
instantaneous rate that corresponds to the gas distribution. 
Moreover, as the right panel of Figure [5] shows, the num- 
ber of star-forming particles in a given simulated galaxy is 
~ 10 2 times larger than that of young stellar particles. This 
ratio could be understood by noting that the observed star 
formation law implies a gas consumption time scale in the 
ISM that is ~ 10 2 times longer than the life times of mas- 
sive stars. Hence, if star formation is implemented by the 
stochastic conversion of gas particles into stellar particles, 
as is the case here, the number of star- forming (i.e., ISM) 
particles is expected to be ~ 10 2 times larger than the num- 
ber of young stellar particles. This ratio will be somewhat 
smaller in starbursts or if gas particles are allowed to spawn 
multiple stellar particles, but under realistic conditions, star 
formation will be sampled substantially better by gas parti- 
cles than by young stellar particles. 

As a result of the above mentioned sampling effects, us- 
ing stellar particles as ionizing sources (as was done in all 
previous work, e.g. , iNagamine et ail |2010| ; iFumagalli et al.l 
l201ll ; lYaiima et al.l I2012T I would underestimate the impact 
of local stellar radiation on the HI distribution for a large 
fraction of simulated galaxies. We therefore use star-forming 
SPH particles instead of young stellar particles as local 
sources of radiation. 

Figure[S]illustrates the difference between the photoion- 
ization rate profiles produced by star-forming SPH particles 
(blue solid curves) and young stellar particles (red dotted 
curves) at z — 3. The left panel of Figure HJ] shows this for 
a halo with M200 = 2 x 10 11 Mq (see the top panels of 
Figure [2}, while the right panel illustrates the results for 
haloes with 10 < log 10 M200 < 10.5 Mg. For this compari- 
son we imposed the same total number of emitted photons 
in the simulation box in both cases. This was done by set- 
ting the photon production rates of individual gas particles 
proportional to their star formation rates (see equation [2]) 
and setting the photon production rate of individual young 
stellar particles proportional to their masses. 

For the massive halo, simulating local stellar radiation 
using SPH particles results in photoionization rate profiles 
similar to those produced by using stellar particles. This 
similarity is mainly due to the agreement between the total 
star formation rate and the rate by which gas is converted 
into young stellar particles (see the left panel of Figure [5]). 
However, this galaxy has ~ 50 times fewer young stellar par- 
ticles than star-forming SPH particles (see the right panel 
of Figure O and as the shaded areas in the left panel of Fig- 
ure [6] show, the scatter in the photoionization rate profile is 
much larger when stellar particles are used. 

For haloes with slightly lower masses, shown in the right 
panel of Figure [(J the rates by which gas is converted into 
young stellar particles are similar to the total star formation 



Table 1. The comoving mean free path of hydrogen Lyman- 
limit photons, A m j p , at different redshifts. From left to right, 
columns respectively show redshift, A m f p in comoving Mpc 
(cMpc) and the references from which the mean fr ee path 
values are taken, i.e. B07: Bolton fc Haehnerd d2007h. F G08: 
IFaucher-Giguere et all ll2008al) , P09: iProchaska et alj [|2009h and 
SC10: ISongaila fc Cowij fcOlOri . 



Redshift 


A m fp 


Reference 




(cMpc) 




2 = 2 


909 ± 252 


FG08, SC10 


2 = 3 


337 ± 170 


FG08, SC10 


2 = 4 


170 ± 15.5 


FG08, P09 


2 = 5 


83.4 ±21.6 


B07, SC10 



rates (see the left panel of Figure [5]). Despite this agreement, 
the median photoionization rates for these two cases differ 
dramatically, being much higher when star-forming gas par- 
ticles are used as sources. We conclude that using stellar 
particles as sources results in the underestimation of the 
ionization impact of local stellar radiation on most of the 
galaxies in the simulation. The intensity of the UVB pro- 
duced by the simulation that uses young stellar particles 
as sources does not agree with the observed UVB intensity, 
while using star-forming particles resolves this issue, after 
correcting for the box size (see £14. 3[) . 

Although in this section we have shown that using 
star-forming particles as sources of ionizing radiation helps 
to resolve the above mentioned sampling issues in post- 
processing RT simulations, we note that doing the same may 
not be a good solution for simulations with sufficient resolu- 
tion to model the cold, interstellar gas phase. Such simula- 
tions can capture the effects of the relative motions of stars 
and gas, e.g., the effect of stars moving out of, or destroying, 
their parent molecular clouds. 



4.3 Stellar ionizing radiation, its escape fraction 
and the buildup of the UVB 

Observations show that the luminosity fun ction of bright 
quasa rs drops sharply at redshifts z > 2 fe.g.. lHopkins et ahl 
l2007h and models for the cosmic ionizing background indi- 
cate that star-forming galaxies dominate th e production of 
hydr o gen ionizing photons at z > 3 (e.g.. lHaehnelt et al.l 
120011 ; Bolton et afll2005t IFaucher-Giguere et al.ll2008al ). In 
the simulations, it should therefore be possible to build up 
the background radiation from the radiation produced by 
star-forming galaxies at z > 3. In this section, we quantify 
the contribution of the ionizing photons that are produced 
in stars to the build-up of the UVB radiation in our simula- 
tions. In addition, we calculate the implied average escape 
fraction that is required to generate the observed UVB. 



4-3.1 Generating the UVB 

Our simulation box is smaller than the typical mean free 
path of ionizing photons, A m f P (see Table [1]). Therefore, the 
IGM gas can receive stellar ionizing radiation from a re- 
gion that is larger than the simulation box. We note that 
using periodic boundaries to propagate the stellar ionizing 
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Figure 7. Left: The simulated UVB photoionization rate produced by stellar radiation in our reference simulation is shown with filled 
circles for different redshifts. The observed mean free paths of ionizing photons have been used to correct for the s mall size of the simulation 
box (see text). The error bars represent the 1 — a errors in these mean free paths. The dashed curve shows the lHaa rdt & Madau (2001) 
UVB photoionization r ates that have been used i n our simulations as the UVB. The observational measurement of the UVB from the 
Lycr effective opacity bv lBolton fc Haehnelj [1200711 is shown using orange diamonds. Right: the average escape fraction of stellar ionizing 
radiation into the IGM calculated based on equation 1131 Our simulation reproduces the observed UVB photoionization rates between 
z = 2 and 5. The implied average escape fractions are 10 — 2 < / CS c < 10~ x between 2 = 2 and 5 and they decrease with decreasing 
redshift below z = 4 



photons is not a good solution to resolve this issue. In addi- 
tion to increasing the computational expense, using periodic 
boundaries for RT would require us to account for the cos- 
mological redshifting of ionizing photons. Moreover, because 
of the small size of our box, stellar photons traveling along 
paths that are nearly parallel to a side of the box may never 
intersect an optically thick absorber, if periodic boundaries 
are being used for RT. However, we can correct for the 
small size of the simulation box by requiring consistency be- 
tween the simulated UV intensities that local sources pro- 
duce in the IGM and existing observations/models of the 
UVB. In order to do that, we assume an isotropic and ho- 
mogeneous universe that is in photoionization equilibrium, 
and A m f p /(l + z) -C c/H(z), so that we can ignore evolution 
and redshifting during the travel time of the photons. Then, 
the volume- weighted mean ionizing flux, from stars is 
given by: 



F* = 



nfp dr — u t X 



mfp , 



(8) 



where u* is the photon production rate per unit comoving 
volume. Moreover, one can express the mean ionizing flux in 
equation Q in terms of the volume-weighted mean ionizing 
flux produced by source that are inside the simulation box, 



nfp cJ r) 



(9) 



where a < 1 is a geometrical factor. In the absence of any 
absorption and for uniform and isotropic distributions of gas 
and sources, a is set by the average distance between two 
random points inside a cube. For a cu be with unit length this 
would yield a = 0.66 (|Robbindll978h . Based on the average 



distance between the low-density gas (e.g., SPH particles 
with n H < 10 -5 cm" 3 ) and sources (i.e., star-forming parti- 
cles) in our simulations, we find that the average value of a 
is 0.66 which varies mildly with redshift from ao = 0.54 at 
z — to as = 0.79 at z = 5. 

As Table [I] shows, our simulation box is much smaller 
than the mean free path of ionizing photons. Therefore, we 
can use I/box <C A m f p in equation © and get: 



(10) 



Using equation (JSJ and (|10|1 , the total stellar UVB radiation 
flux can thus be written as: 

A m f P 



Ft, 



a L b 



(11) 



Therefore, the volume-weighted photoionization rate due to 
radiation produced by stars in the simulation, T™, which is 
close to the median photoionization rate in low-density gas, 
can be used to calculate the implied UVB photoionization 
rate after correcting for the small box size of the simulations: 

„ A m fp 



a Lb 



(12) 



As mentioned earlier, our simulations show that the 
photoionization rate from local stellar radiation approaches 
a density independent rate at low densities. We take 
this photoionization rate as F 1 * 11 . In addition, we use 
a compilation of available Lyman-limit mean free path 
measurem ents at different redshif t s from the literature 
(i.e., f r om Bolton fc HaehneltJ 2007 : Faucher-Giguere et al.l 
l2008al : IProchaska et all 120091 : lsongaila fc Cowid l20ld : see 
Table QJ. After converting the Lyman-limit mean free paths 
into the typical mean free path of ionizing photons with our 
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assumed stellar spectrurr0, we use equation (|12[) to derive 
the implied UVB photoionization rate. 

Figure [7] shows the predicted contribution of stellar ra- 
diation to the UVB (filled circles). The error bars reflect the 
quoted error in the mean free path measurements. For com- 
parison, the observational measurement of the UVB from 
the Lya effe ctive opacity bvlBolton fc Haehnehj (|2007t ) and 
the modeled lHaardt fc Madau 1 200 if ) UVB photoionization 
rates are also shown by orange diamonds and green dashed 
curve, respectively. Both the observational measurements 
and modeled UVB intensities are in good agreement with 
our simulation results for z > 2. However, their UVB inten- 
sity is slightly higher than ours at z = 2. The reason for 
this could be the absence of radiation from quasars in our 
simulations. 



HI absorbers with different strengt hs in the spec- 
tra of background quasars (e.g., Kim et al. 20021; 

[ Peroux et~ai]|2005l; Ib'Meara et al.ll2007l; iNoterdaeme et all 
20091; IProchaska et all 120091; IProchaska fc Wolfe! 12009^ 
lo'Meara et all |2012| ; fNoterdaeme et all l2012t ). The HI 
column density distribution function (CDDF) is defined as 
the number of systems at a given column density, per unit 
column density, per unit absorption length, dX: 



f(N m ,z) 



d 2 n 



d 2 n H(z) 



dN m dX dN m dz H (1 + z) 2 



(14) 



In order to study the effect of local stellar radiation on 
the HI CDDF, we project our simulation boxes on a two- 
dimension al grid and use this to calculate the column den- 
sities (see lRahmati et al]|2012i for more details). 



4-3.2 Average escape fractions 

The average star formation rate density in our simulations is 
in go od agreement with the observed cosmic star formation 
rate ijSchaveet alj|20ich . Therefore, the good agreement be- 
tween the UVB intensities in our simulation and the ones 
inferred from the observations suggests that at z > 3 the 
average escape fractions of stellar ionizing photons in our 
simulations are also reasonable. However, as we will discuss 
in Appendix IB21 the structure of the ISM is unresolved in 
our simulations. Therefore, the fact that the produced es- 
cape fractions are reasonable may be coincidental. 

Since in our simulations both the intensity of stel- 
lar radiation in the IGM and the photon production rate 
are known, we can measure the mean star formation rate- 
weighted escape fraction of ionizing photons from galaxies 
into the IGM (see Appendix [C]): 



10" 



10- 



2.9 x 10- 18 cm 2 
1 



0.15 M yr^cMpc" 3 

JM _1 / £box V 1 / 1 
0.7/ VlOcMpc/ V 



+ z 



(13) 



The implied escape fractions are shown in the right panel 
of Figure [7] The simulated escape fractions are 10" 2 < 
/esc < 10 _1 between z = 2 and 5 and they decrease 
with decreasing redshift below z = 4. This result is con- 
siste nt with previous observational and theoretical studies 
(e.g.. IShaplev et al.ll2006l; ISchavell20of3 ; iGnedin et af1l2008l ; 
iKuhlen fc Faucher-Giguerdl2012i ). 



4-4-1 Local stellar radiation and the HI CDDF at z = 3 

The simulated HI CDDF at z = 3 is shown in the left panel 
of Figure[8]for different ionizing sources. The blue solid curve 
shows the HI CDDF in our reference simulation which in- 
cludes local stellar radiation, the UVB and recombination 
radiation. For comparison, the HI CDDF without local stel- 
lar radiation (i.e., only including the UVB and recombina- 
tion radiation) is shown with the green dashed curve. As the 
ratio between these two HI CDDFs in the top section of the 
left panel in Figurc[8]illustrates, the effect of local stellar ra- 
diation increases with HI column density and reaches a ~ 0.5 
dex reduction at A'hi > 10 21 cm -2 . For HI column densities 
lower than 10 17 cm" 2 on the other hand, the HI CDDF is 
insensitive to local stellar radiation. These trends are con- 
sistent with previ ous analytic arguments dMiralda-Escudej 
20051; Schavell2006l) and numerical simulations performed by 



Fumagalli et all l|201lT l 



As we discussed in £14.21 using star-forming SPH par- 
ticles as ionizing sources (i.e., our reference model) results 
in a better sampling than using young stellar particles as 
sources. We showed that if one uses young stellar particles 
as ionizing sources, the impact of local stellar radiation will 
be under-estimated in a large fraction of low-mass haloes in 
our simulations (see Figure [5] and [5]). To illustrate this differ- 
ence, the HI CDDF for the simulation in which young stellar 
particles are used, is shown with the red dotted curve in the 
left panel of Figure [S] Indeed, the impact of local sources 
on the HI CDDF is much weaker if we use young stellar 
particles as ionizing sources. This may partly explain why 
lYaiima ~et al. (2012) found that local sources did not affect 
the HI CDDF significantly. 



4.4 The impact of local stellar radiation on the 
HI column density distribution 

The observed distribution of neutral hydrogen is 
often quantified by measuring the distribution of 



2 Because the effective hydrogen ionization cross section of stellar 
ionizing photons that we use in this work, cr* = 2.9 X 10 — 18 cm 2 , 
their typical mean free paths are ~ 2 times longer than the mean 
free path of Lyman-limit photons which have hydrogen ionizing 
cross sections <tq = 6.8 X 10~ 18 cm 2 . 



4-4-2 The impact of the unresolved ISM 

In our simulations, the ISM is modeled by enforcing a poly- 
tropic equation of state on SPH particles with densities 
ro H > 10 _1 cm~ 3 . This means that our simulations do not 
include a cold (T <C 10 4 K) interstellar gas phase. This sim- 
ple ISM modeling introduces uncertainties in the hydrogen 
neutral fraction calculations of the dense regions in our simu- 
lations. In addition, for gas with n H 3> 10 -2 cm -3 the mean 
free path of ionizing photons is unresolved (see Appendix 

To estimate the impact of the structure of the ISM on 
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log 10 ( N HI [cm" 8 ] ) log 10 ( N HI [cm" 8 ] ) 

Figure 8. Left: The HI column density distribution function (CDDF) at z = 3 with different ionizing sources. The blue solid curve 
shows our reference simulation which includes the UVB, local stellar radiation (LSR) and recombination radiation and the green dashed 
curve indicates the simulation without local stellar radiation (i.e., with the UVB and recombination radiation). While in the reference 
simulation star-forming SPH particles are used as ionizing sources, the HI CDDF that is shown with the red dotted curve indicates a 
simulation in which young stellar particles are used as sources. Using star-forming SPH particles as sources lowers the HI CDDF by 
~ 0.5 dex for Nm > 10 21 cm -2 . However, using young stellar particles, which results in sampling issues (see ^4.21 . has a weak impact 
on the HI CDDF. For calculating the HI CDDF, the neutral gas is assumed to be fully atomic (i.e., no H2). Right: The HI CDDF at 
z = 3 in the presence of local stellar radiation (i.e., star-forming SPH particles as sources) for different assumptions about the ISM. 
The blue solid curve shows our reference simulation in which all H atoms contribute to the absorption but molecular hydrogen does not 
contribute to the HI CDDF. The orchid dot-dashed curve shows a porous ISM model (see the text) where molecular hydrogen does not 
absorb ionizing radiation during the RT calculation (it is assumed to have a very small covering fraction). In order to reproduce the 
observed UVB intensity in this model, the local stellar radiation has been reduced by a factor of 3 compared to the fiducial model. The 
red long-dashed curve is identical to the orchid dot-dashed curve (porous ISM) but molecular hydrogen is assumed to dissociate into 
atomic hydrogen before calculating the HI CDDF. At higher HI column densities (TVhi > 10 21 cm~ 2 ) the HI CDDF is highly sensitive 
to the assumptions about the unresolved ISM. 



the HI CDDF and the effect of local sources, we introduce a 
porous ISM model in which we assume that molecular hydro- 
gen is confined to clouds with such small covering fractions 
that we can ignore the m when performing the RT. Following 
iRahmati etafl j|2012r). we use the observa tionally inferred 
pressure law of lBlitz fc Rosolowskvl (|2006l ) to compute the 
molecular fraction (see Appendix [A} , but now we not only 
subtract the H2 fraction when projecting the HI distribution 
to compute the HI CDDF, we also subtract it before doing 
the RT calculation (note that in the fiducial case we assumed 
H2 fractions to be zero during the RT calculation) . The con- 
version of diffuse atomic hydrogen into compact molecular 
clouds which do not absorb ionizing photons should facili- 
tate the propagation of photons from star-forming regions 
into the IGM. Indeed, we find that replacing our reference 
uniform ISM model with the porous ISM model increases 
the resulting UVB photoionization rate by a factor of ~ 3. 
Therefore, we decreased the photon production rate of local 
stellar radiation by a factor of 3, such that the model with 
a porous ISM yields the same UVB intensity that is gener- 
ated by our reference simulation and is in agreement with 
the observed UVB. This factor of 3 reduction could be in- 
terpreted as reflecting absorptions in molecular clouds that 



are hosting young star^f]. One should note that despite the 
enforced agreement between the generated UVBs, the stel- 
lar photoionization rates at high and intermediate densities 
are different in the two simulations. 

The right panel of Figure [8] shows that for A^hi < 
10 21 cm -2 the HI CDDF in the simulation with a porous 
ISM (orchid dot-dashed curve) is almost identical to the HI 
CDDF in our reference simulation (blue solid curve). This 
suggests that for these column densities the predicted im- 
pact of local stellar radiation on the HI CDDF is robust to 
uncertainties regarding the small scale structure of the ISM. 
Instead, its impact is controlled by the amount of radiation 
that is propagating through the properly resolved interme- 
diate densities towards the IGM, which is constrained by the 
observed UVB intensity. Note, however, that the impact of 
local sources on the high HI column density part of the dis- 
tribution (i.e., JVhi > 10 21 cm -2 ) is in fact highly dependent 



3 Note that this is not equivalent to / csc = 1/3. The radiation 
that leaves star-forming regions is still subject to significant ab- 
sorption before it reaches the IGM. 
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Figure 9. The ratio between the HI CDDF with and without 
local stellar radiation at different redshifts. For calculating the 
HI CDDF, the neutral gas is assumed to be fully atomic (i.e., 
no H2). The impact of local stellar radiation decreases the HI 
CDDF by up to 1 dex. The impact of local stellar radiation on 
LLSs increases with redshift. 

on the assumptions that are made about the physics of the 
ISM. 

It is also possible that the intense stellar ionizing ra- 
diation within the ISM will effectively increase the HI col- 
umn densities by dissociating hydrogen molecules. However, 
we implicitly neglected this effect in our simulation with a 
porous ISM by assuming that molecular clouds are not af- 
fected by local stellar radiation. To put an upper limit on 
the impact of this effect, one can assume that all molecular 
clouds are completely dissociated by the absorption of stel- 
lar radiation. The result of this exercise is shown by the red 
long-dashed curve in Figure [8] This shows that accounting 
for H2 dissociation could reduce (or even reverse) the im- 
pact of local stellar radiation only at the very high column 
density end of the HI CDDF (i.e., Nm > 10 21 cm" 2 ). 

We stress that none of our ISM models are realistic. 
However, by considering very different models, we can nev- 
ertheless get an idea of the possible impact of our simplified 
treatment of the ISM and we conclude that our results are 
relatively robust for Am <C 10 21 cm -2 . 

4-4-3 Evolution 

The evolution of the impact of local stellar radiation on the 
HI CDDF is illustrated in Figure [9] for our fiducial model. 
Each curve in this figure shows the ratio between the HI 
CDDF with and without local stellar radiation at a given 
redshift. To avoid the uncertainties about the conversion of 
atomic gas into H2, the HI CDDFs are computed assuming 
that the neutral gas is fully atomic (i.e., no H2). For all red- 
shifts local stellar radiation has only a very small impact 
on the HI CDDF for N m < 10 17 cm" 2 but significantly re- 
duces the abundance of systems with Am > 10 21 cm -2 . The 
impact of local stellar radiation increases with redshift for 



10 18 < iVm < 10 21 cm 2 . While local sources significantly 
reduce the HI CDDF for A m > 10 17 cm" 2 at z = 5, their ef- 
fect only becomes significant for Am > 10 21 cm" 2 at z = 0. 
This might be attributed to decrease in the proper sizes of 
galaxie s with redshift. 

In iRahmati et alJ l|2012h we used simulations that in- 
clude only the UVB and recombination radiation to show 
that for 10 18 < N m < 10 20 cm" 2 , the HI CDDF does not 
evolve at z < 3 and increases with increasing redshift for 
z > 3. Figure [10] shows that if we also include local stellar 
radiation, the weak evolution of the HI CDDF in the Lyman 
Limit range extends to even higher redshifts (i.e., z < 4). 

The predicted HI CDDFs are compared to observations 
in Figure \W\ We note that the HI CDDF predicted by our 
reference simulation is not converged with respect to the 
size o f the simulation box (see Appendix B in lRahmati et ahl 
I2012T ). Since the photoionization caused by local stellar ra- 
diation only exceeds the UVB photoionization rate close to 
galaxies (see Figure [4]), we can assume that the impact of 
local stellar radiation on the HI CDDF is independent of 
the size of the simulation box. Therefore, for all the curves 
shown in Figure [10] we have corrected the box size effect 
by multiplying the CDDF of the fiducial L = 6.25/i -1 Mpc 
box by the ratio of the CDDF in a converged simulation 
(with L — 50ft" *Mpc) and in our reference simulations with 
L = 6.25ft _1 Mpc, both in the absence of local stellar radia- 
tion (i.e., with the U VB and recombination radiation). 

As we showed in I Rahmati et all ([2012), the simulated 
HI CDDFs in the presence of the UVB and recombination 
radi ation is in a reason able agreement with observations (see 
also lAltav et alJl201ll ) with a small deviation of ~ 0.2 dex. 
But as the top section of Figure flOl illustrates, the addition 
of local stellar radiation increases this deviation to > 0.5 
dex. While the difference between the predicted HI CDDFs 
and observations is larger at the highest HI column densities 
(i.e., Am > 10 21 cm" 2 ), we have seen that in this regime the 
effect of local sources is sensitive to the complex physics of 
the ISM, which our simulations do not capture. These very 
high HI column densities are also sensitive to the strength 
and the details of different feedback mechanisms (see Altay 
et al. in prep.). The small but significant discrepancy in the 
LL and weak DLA regime is therefore more interesting. It 
is important to confirm this discrepancy with larger simu- 
lations, so that a correction for box size will no longer be 
necessary, but this requires more computing power than is 
presently available to us. 



5 DISCUSSION AND CONCLUSIONS 

The column density distribution function (CDDF) of neutral 
hydrogen inferred from observations of quasar absorption 
lines is the most accurately determined observable of the 
distribution of gas in the high-redshift Universe. Moreover, 
the high column density absorbers (Nm > 10 17 cm" 2 ) arise 
in gas that is either a lready part of the ISM o r will soon 
accrete onto a galaxy (|van de Voort et al . 2012) and hence 
these systems directly probe the fuel for star formation. 

To predict the distribution of high column density ab- 
sorbers, it is necessary to combine cosmological hydrody- 
namical simulations with accurate radiative transfer (RT) 
of ionizing radiation. Because of the cost and complexity 
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Figure 10. The HI CDDF at different redshifts. Curves show 
predictions for the reference simulation in the presence of the 
UVB, diffuse recombination radiation and local stellar radiation 
after correcting for the box size (see the text). For calculating 
the HI CDDF, the neutral gas is assumed to be fully atomic (i.e., 
no H2). The observational data points represent a compilation of 
various quasar absorption line observa tions at high redshifts (i.e., 
2 = [1.7,5.51) t a ken fr om lPeroux et al.l 1120051 1 with 2 = [1.8,3.51, 
lO'Meara et al.l J2007tl with z = [1.7, 4.51. iNoterdaeme et all 
teOOSft withz = [2.2, 5.5] and iProchaska fc Wolfe! l l2009t) with 
z = [2.2, 5.5]. The orange filled circles s how the best - fit bas ed on 
the low-redshift 21 cm observations of IZwaan et al. I ||2005| 1. The 
top section shows the ratio between the HI CDDFs at different 
redshifts and the predicted CDDF at 2 = 3. While simulations 
without local stellar radiation are in reasonable agreement with 
the o bserved HI CDDF (not shown here, but see lRahmati et al.l 
2012), in the presence of local stellar radiation the simulated HI 
CDDFs of LLSs deviate from observations by a factor of fa 2. 
Local stellar radiation weakens the evolution of the HI CDDFs at 
10 17 < N m < 10 20 cm" 2 . 



associated with RT calculations that include many sources, 
nearly all studies have only considered ionization by the ul- 
traviolet background radiation (UVB), whose intensity can 
be inferred from observations of the Lya forest. However, 
analytic arguments suggest that the radiation field to which 
high column density absorbers are exposed is typically dom- 
inate d by local stellar sources (jMiralda-Escud 6 2005; Schaye 
2006) . It is therefore important to investigate whether the re- 



markable success of simulations that consider only the UVB, 
such as the agreement with the observed CDD F over 10 or- 
ders of magnitude in column density at z = 3 (lAltav et al 
2011T ) as well as its evolution down to 2 = (jRahmati et al 



20121 ) . is compromised when local sources are included. Here 



we add ressed this question by repeating some of the simula- 
tions oflRahmati et all &20m with our RT code T RAPHIC 
jPawlik fc Schavd200SU201ll ; lRahmati et al.ll2012l ). but this 
time including not only the UVB and diffuse recombination 
radiation, but also local stellar sources. 

In agreem ent with the analytic predictions of 

iMiralda-Escudi l|2005l ) and ISchavd (|2006l ). we found that 
local stellar radiation is unimportant for JVhi <C 10 17 cm" 2 
and dominates over the UVB for high column density ab- 
sorbers. For all redshifts considered here (i.e., z < 5), lo- 
cal sources strongly reduce the abundance of systems with 
A^hi > 10 21 cm" 2 . The impact of local sources increases with 
redshift for 10 18 < N m < 10 21 cm" 2 . At 2 = 5 the CDDF 
is substantially reduced for TVhi ^> 10 17 cm" 2 , but at z = 
the effect only becomes significant for Nm > 10 21 cm -2 . 
As a result, the r emarkable lack of evo lution in the CDDF 
that we found in iRahmati et al. (2012) for 2 = — 3, and 
which is also observed, extends to z = 4 if local sources are 
taken into account. On the other hand, the agreement with 
the observed 2 ~ 3 CDDF is not quite as good as before, 
with the simulations underpredicting the rates of incidence 
of 10 19 < N m < 10 21 cm absorbers by factors of a few. 
However, because of the large corrections that we had to 
make because of the small size of the simulation box used 
to study the effect of local sources (6.25/i _1 Mpc), this dis- 
crepancy will have to be confirmed with larger simulations. 
Moreover, we did not account for possible hydrodynamical 
effects that might be caused by extra heating due to ionizing 
radiation from local sources. This process might change the 
distribution of gas that is affected by local stellar radiation 
and requires further investigation. 

We found that the average photoionization rate due 
to young stars in high-density gas is weakly dependent on 
the gas density and is ~ 10 _13 s _1 . We showed analytically 
that this rate follows directly from the imposed (and ob- 
served) Kennicutt-Schmidt star formation law if we assume 
that most of the ionizing photons that are produced by star- 
forming gas are absorbed on scales < kpc. However, in reality 
we expect the photoionization rate in the ISM to fluctuate 
more strongly than predicted by simulations like ours, which 
lack the resolution required to model the cold, interstellar 
phase. 

Indeed, the spatial resolution that is required for ac- 
curate RT of ionizing radiation through the ISM is several 
orders of magnitude higher than the smallest scales accessi- 
ble in current cosmological simulations. This makes tackling 
this problem in the near future hardly feasible and poses 
a difficult challenge for studying the impact of local stellar 
radiation on the distribution of HI in and around galaxies. 
Fortunately, one can circumvent part of the problem by tun- 
ing the production rate of ionizing photons (which is equiv- 
alent to adjusting the escape fraction of ionizing photons 
from the unresolved ISM) such that the models reproduce 
the observed mean photoionization rate in the IGM (after 
subtracting the contribution from quasars). In other words, 
if one knows the amount of ionizing radiation that is required 
to reach the IGM, its ionization impact on the intervening 
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gas can be determined even if we cannot predict what frac- 
tion escapes from the immediate vicinity of the young stars. 

We adopted this approach but found that tuning was 
unnecessary for our reference simulation at z ~ 3. More- 
over, since our simulations also yield star form ation histories 
that are in good agreement with observations jSchave et al.l 
|2010T ). we used them to constrain the implied star forma- 
tion rate weighted mean escape fraction that relates the 
predicted star formation rate density to the intensity of the 
UVB. We found that the average escape fraction in our sim- 
ulations is 10~ 2 — 1CP 1 at z — 2 — 5, which agrees with pre- 
vious constramtson_the escape fraction from obser vations 
(e.g.,ls"naplev et alj|2006t ) and theoretical work (e.g.,|Scha 
l200d : lGnedin et aUbooa iKuhlen fc Faucher-Gigu^rdl20l" 



The limited spatial resolution of cosmological simula- 
tions mandates the use of simplified models for the structure 
of the ISM. To estimate the impact of such subgrid models 
on the CDDF and on the effect of local sources, we var- 
ied some of the underlying assumptions. In particular, we 
considered a porous ISM model which assumes that molec- 
ular hydrogen is confined to clouds with such small covering 
fractions that we can ignore them when performing the RT. 
We also considered a model which assumes that all molecu- 
lar clouds are completely dissociated, but not ionised, by the 
absorption of stellar radiation. Although none of these mod- 
els are realistic, we used them to estimate the potential im- 
pact of our simplified treatment of the ISM. We found that 
provided that we rescale the source luminosities so that the 
different models all reproduce the observed background radi- 
ation, the models predict nearly the same HI CDDFs in the 
regime where the absorbers are well resolved. We therefore 
concluded that our results on the effect of local sources are 
relatively robust for TVhi <C 10 21 cm -2 , but that their pre- 
dicted impact is highly sensitive to the assumptions about 
the ISM for N m > 10 21 cm" 2 . 

Different studies have found qualitatively different re- 
sults for the i mpact of local st e llar ra diation on the CDDF. 
For instance, iFumagalli et al.l (|201 lh used relatively high- 
resolution zoomed simulations of individual objects to 
demonstrate that local stellar radiation significantly reduces 
the abundance of high column density absorbers. On the 
other hand, some studies using cosmological simulations 
similar to those presented here (with roughly the same res- 
olution) found that lo cal stellar radiation h a s a negligible 
impact on the CDDF l|Nagamine et al.ll2OI0l : lYai ima et al.l 
l2012h . 

We found that difference in the resolutions of the simu- 
lations that were used in these previous studies may explain 
their inconsis tent findings. Using star particles as source s, as 
was done by l|Nagamine et al.ll2010l ; lYaiima et al.ll2012T ). we 
also found that local sources have a negligible impact on the 
abundance of strong HI systems. However, we demonstrated 
that it is possible to dramatically improve the sampling of 
the distribution of ionizing sources by using star-forming gas 
particles (i.e., gas with densities at least as high as those 
typical of the warm ISM), thus effectively increasing the 
resolution without modifying the time-averaged production 
rate of ionizing photons. We adopted this strategy in our 
fiducial models and found, as summarized above, that the 
radiation from local sources significantly affects the high col- 
umn densitv_end_ofth«_CDDR result is in agreement 
with IFumagalli et all (|201ll ). lGned in (2010) and analytic es- 



timates of iMiralda-Escudel (|2005l ) and ISchavej (|2006h . and 
confirms that poor sampling of the distribution of ionizing 
sources can lead to an under-estimation of the impact of 
local stellar radiation. 

Further progress will require higher resolution simula- 
tions and, most importantly, more realistic models for the 
ISM. In the near future it will remain unfeasible to accom- 
plish this in a cosmologically representative volume. Until 
this challenge is met, predictions for the escape fractions of 
ionizing radiation averaged over galaxy populations should 
be considered highly approximate. Predictions for the abun- 
dances of LL and weak DLA systems based on models that 
neglect local sources of stellar radiation should be inter- 
preted with care, particularly for z > 2. Predictions for the 
CDDF in the strong DLA regime (Nm > 10 21 cm" 2 ) must 
be considered highly approximate at all redshifts. 
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APPENDIX A: HYDROGEN MOLECULAR 
FRACTION 

To account for the effect of molecu lar hydrogen, we adopt 
the same H2 conversion relation that lAltav et alj (|201lt ) used 
to suc cessfully reproduce the HI CDDF high-end cut-off. We 
follow [Blitz fc Rosolowskvl (120061 ) and adopt an observation- 
ally inferred scaling relation between the gas pressure and 
the ratio between molecular and total hydrogen surface den- 
sities: 



Pcxt 



(Al) 



where J? mo i = Eh 2 /£hi, -Pext is the galactic mid- plane pres- 
sure, a = 0.92 and Po/kb = 3.5 x 10 4 cm -3 K. Furthermore, 
if we assume J? mo i gives also the local mass ratio between 
molecular and atomic hydrogen, the fraction of gas mass 
which is in molecular form can be written as 
Mh 2 M h , 1 



/h 2 = 



Mtot Mh 2 + Mm 1 + 7?" 



(A2) 



The last equation also assumes that at high densities ioniza- 
tion is not dominant and the gas is either molecular or neu- 
tral. In our simulations we model the multiphase ISM by im- 
posing an effective equation of state with pressure P oc p 7oft 
for densities hh > n^, where n£j = 0.1 cm -3 which is nor- 
malized to P*/k b = 1.08 x 10 3 cm -3 K at the threshold. 
Therefore we have: 



pA^ 

n„ 



(A3) 



To make the Jeans mass and the ratio of the Jeans length to 
the SPH kernel independent of the density, we use y e a =4/3 
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l|Schave fc Dalla Vec chia 2008). Consequently, after combin- 
ing equations (|A1|I . (|A2f) & (|A3f) the fraction of gas mass 
which is converted into H2, /h 2 can be written as: 

= + * (if)"') ^ 
where A = (P*/P )~ a = 24.54 and /3 = a j cB = 1.23. 



APPENDIX B: RESOLUTION EFFECTS 

Bl Limited spatial resolution at high densities 

At high densities, ionizing photons are typically absorbed 
within a short distance from the sources. The propaga- 
tion of ionizing photons in these regions is therefore con- 
trolled by distribution of gas on scales that are comparable 
to the very short mean free path of ionizing photons. In 
Figure IB11 the mean free path of stellar ionizing photons, 
Amfp = l/( n m <^*fl is plotted as a function of density for 
our reference simulation at * = 3. The blue dashed curve in 
Figure [Bl] shows the result when only the UVB and recom- 
bination radiation are included. The effect of adding local 
stellar radiation is shown by the red dot-dashed curve. For 
comparison, also the result for a constant UVB radiation 
(i.e., optically thin gas) is illustrated with the green dotted 
curve. The ionization by local stellar radiation decreases the 
hydrogen neutral fraction at n H > 10 -2 cm -3 which in turn 
increases the mean free path of ionizing photons making it 
comparable to the optically thin case. 

The simulations cannot provide any reliable information 
about the spatial distribution of gas and ionizing sources on 
scales smaller than the resolution limit (i.e., the typical dis- 
tance between SPH particles) which is shown by the black 
dotted line in Figure IB1I At densities for which the mean 
free path of ionizing photons is shorter than the spatial reso- 
lution, the RT results may not be accurate since all the pho- 
tons that are emitted by sources are absorbed by their imme- 
diate neighbors. This happens at densities n H > 10~ 2 cm -3 
without local stellar radiation (blue dashed curve in Figure 
IB1|) and at densities n H > 10 _1 cm -3 with local stellar ra- 
diation (red dot-dashed curve in Figure [Bljl . On the other 
hand, Rt is irrelevant if the gas is highly neutral, which is 
predicted to be the case at slightly higher densities (see Fig- 
ure E]). 



B2 The impact of a higher resolution on the RT 

In iRahmati et al l (|2012h we showed that the self-shielding 
limit is not very sensitive to the resolution of the simulation. 
On the other hand, the small scale structure of the ISM may 
significantly change the propagation of stellar radiation and 
their impact on the HI distribution. While this suggests that 
the impact of local stellar radiation might be sensitive to 
the resolution of the underlying simulation, our approach of 



4 (7 -it is the spectrally averaged hydrogen photoionization cross 
section for stellar ionizing photons. For the spectral shape of the 
stellar radiation, we adopt a blackbody spectrum with T^b = 

5 X 10 4 K. 
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Figure Bl. The mean free path of stellar ionizing photons as a 
function of hydrogen number density for the reference simulation 
at z = 3. The red dot-dashed curve shows the result when all 
radiation sources, i.e., UVB, local stellar radiation (LSR) and 
recombination radiation (RR), are included and the blue dashed 
curve shows the result when local sources of radiation are not 
present. The green dotted curve shows the mean free path in 
a simulation that assumes a completely optically thin gas. The 
typical distance between SPH particles as a function of density 
is shown by the black dotted line. The Jeans length for a given 
density is illustrated by the black dashed line. The colored lines 
indicate the medians and shaded regions around them show the 
15% - 85% percentiles. 



tuning the escaped stellar radiation such that it can generate 
the desired UVB, circumvent most resolution effects. 

To study the impact of resolution on the propagation of 
stellar ionizing photons in the ISM and their escape to the 
IGM, we performed a simulation similar to our reference 
simulation but using 8 times more dark matter and SPH 
particles (i.e., using 256 3 dark matter particles and the same 
number of SPH particles whose masses are 8 times lower 
than in our reference simulation) . Figure IB2I shows that at 
z = 3, the photoionization rate of stellar radiation in the 
higher resolution simulation (red dashed curve) is qualita- 
tively similar to that obtained from our reference simulation 
(blue solid curve). As expected, the stellar photoionization 
rate at the highest densities is in agreement with our ana- 
lytic estimate in 321 However, because of the shorter inter- 
particle distances in the higher resolution simulation, the 
stellar photoionization rate peaks at slightly higher densi- 
ties. The biggest difference between the stellar photoioniza- 
tion rates of the high-resolution and our reference simulation 
occurs at the lowest densities. Despite having a ~ 2 times 
higher total star formation rate, the high-resolution simula- 
tion results in a UVB which has a ~ 5 times lower photoion- 
ization rate. This means that the effective escape fraction in 
the high-resolution simulation is ~ 10 times lower than for 
the reference simulation which has / csc ~ 10~ 2 . 

Finally, we emphasize that due to our simplified treat- 
ment of the ISM, it is not even clear whether the results 
become more realistic with increasing resolution. 



© 0000 RAS, MNRAS 000, 000-000 



The impact of local stellar radiation on the HI CDDF 



(|C1|) . the escape fraction becomes: 




-5 -4 -3 -2-10 1 
l°g 10 [ n H (cm" 3 )] 

Figure B2. Comparison between the stellar photoionization rate 
profiles in the reference simulation and a simulation with 8 times 
higher resolution at z = 3. The colored lines indicate the medians 
and shaded regions around them shows the 15% — 85% distribu- 
tions. 



APPENDIX C: CALCULATION OF THE 
ESCAPE FRACTION 

Using the equilibrium photon production rate per unit star 
formation rate from equation ([2]), the total production rate 
of ionizing photons in the simulation box, A/" 7 ,*, can be cal- 
culated as a function of the comoving star formation rate 
density within the simulation box, p t , : 



Noting that the typical time the ionizing photons spend in- 
side the simulation box is ~ aLbo*/ (l + z)/c, and letting / csc 
be the star formation rate-weighted mean fraction of ioniz- 
ing photons that escape into the IGM, one can calculate the 
comoving equilibrium photon number density in the IGM: 



£Lx C (1 + Z) 



(C2) 



where we ignored absorptions outside the host galaxies be- 
cause ibex < A mfp . 

The comoving number density of ionizing photons in 
the simulation box is related to the effective hydrogen pho- 
toionization rate they produce in the IGM, r* n : 

" 7 = C <7*(l + Z )3' (C3) 

where ct* is the effective hydrogen ionization cross section for 
stellar photons. Combining equations (|C2[) and (|C3|l yields 
the effective escape fraction of ionizing photons from stars 
into the IGM: 



o * a My,*(l + z) 2 
After putting numbers in equation (IC4|) and using equation 
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